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Abstract 

An analytic model of the shear modulus applicable at temperatures up to melt and 
at all densities is presented. It is based in part on a relation between the melting 
temperature and the shear modulus at melt. Experimental data on argon are shown 
to agree with this relation to within 1%. The model of the shear modulus involves 
seven parameters, all of which can be determined from zero-pressure experimental 
data. We obtain the values of these parameters for 11 elemental solids. Both the ex- 
perimental data on the room-temperature shear modulus of argon to compressions 
of ~ 2.5, and theoretical calculations of the zero-temperature shear modulus of alu- 
minum to compressions of ~ 3.5 are in good agreement with the model. Electronic 
structure calculations of the shear moduli of copper and gold to compressions of 2, 
performed by us, agree with the model to within uncertainties. 
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1 Introduction 

A reliable model of the adiabatic (isentropic) shear modulus, G, of a polycrystalline solid 
at temperatures to T m , the melting temperature, and up to megabar pressures is needed for 
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many applications, including the modeling of plastic deformation at extremes of pressure 
and temperature, numerical calculations of elastic and shock wave propagation, and even 
calculations of the oscillations of low-mass astrophysical objects. 

Adiabatic elastic properties are generally determined by ultrasonic wave-speed mea- 
surements, which are usually made in the low-pressure regime. Zero-pressure experimental 
data have been accumulated on single- crystal elastic constants, together with polycrys- 
talline averages, at temperatures from T = to nearly T m for Ag (to within 60 °K of T m ) 
0, Au (to within 60 °K of T m ) @, Ge (to within 90 °K of T m ) @, and V (to within 80 °K 
of T m ) [§. The data run from T = to T m for Al §, Ar f§, Bi §, Cd [§, Cs §, Cu 
nj, In (TJ, K |g, Na |T|, Nb 0, Ne |T§, Pb Sn §, Ta |T§, Te 0, Xe @, 



and Zn 119 



On the theoretical side, it is possible to calculate singe-crystal elastic constants as 
a function of compression at zero temperature from electronic-structure theory. Such 
calculations were done by Straub et al. [[ZD] for Cu, Christensen et al. P]J for Mo and 

m 



W, Soderlind et al. J22] for Fe, and Soderlind and Moriarty |23) for Ta. With known 
interatomic potentials it is possible to calculate the temperature dependence of the elastic 
constants by computer simulation techniques, as demonstrated by the calculations for Na 
24 1, Mg f2q| , and Cu [26]. Bounds on the shear modulus, G, can be calculated from 
the single-crystal elastic constants for any crystal class [p7 |, and for a cubic crystal the 
polycrystalline shear modulus can be calculated exactly using the Kroner cubic equation 



Guinan and Steinberg j29| modeled the zero-temperature shear modulus as G = Go + 
G' P (p /p) 1/3 , where G' is the pressure derivative of G at zero pressure and p is density. 
This functional form was chosen so that G ~ p 4//3 as p — ► oo, the correct asymptotic 
behavior albeit with a prefactor which does not generally coincide with that given by 
the one-component plasma model for G. Preston and Wallace ]30| proposed a model for 
the temperature dependence of the shear modulus at any density, but left the density 
dependence itself arbitrary. The dependence of the shear modulus on both density and 
temperature has also been discussed by Anderson ]3l] . 



In this paper we develop an analytic model for the density and temperature dependence 
of the shear modulus by combining four key elements. First is a simple but accurate 
relation between the density, the melting temperature as a function of density, T m (p), 
and the shear modulus along the solidus |32|, [33| . Second is the Preston- Wallace model 
for the shear modulus. Third is an analytic model for the Gruneisen parameter [34 that is 
used in conjunction with the fourth ingredient, the Lindemann criterion |35[, to generate 
an analytic expression for T m (p). 



2 A relation between shear modulus and melting 
temperature 

The melting temperature and shear modulus along the solidus approximately satisfy the 
relation 

G(p,T m (p)) = G{p ve{ ,T m (p rc{ )) 

P T m(p) PvcfT m (p rci ) 



where p re f is a reference density. This relation is the foundation of our model for the 
shear modulus, so we provide theoretical justification for it following two approaches: 
the theory of dislocation-mediated melting |3B|, and the theory of a Debye solid (in 
which it derives as a consequence of the proportionality of G to the square of the Debye 
temperature). The relation is shown to agree very well with shear modulus data on argon, 
the only data available for such a comparison. 



2.1 Two derivations of relation (1) 

It follows from our model of melting as a dislocation-mediated phase transition that the 
rolo/tj ion 

1 - v{T m )/2 G(T m )v(T m ) A ( o? \ 
B m 1 - v(T m ) Hz-1) 8*\4Pd(T m ))- [Z) 

holds at any pressure. Here b is the magnitude of the Burgers vector, v is the Poisson 
ratio, v is the Wigner-Seitz volume, A = b 3 /v is a geometric factor characterizing the 
lattice, z is the coordination number, and d(T m ) is the dislocation density at melt. Note 
that the factors A and ln(z — 1) explicitly account for the influence of crystal structure on 
melting. The value of A is 3a/3/4 pa 1.30 for body-centered cubic (bcc), and V2 pa 1.41 



for face-centered cubic (fee) and ideal (c/a = y 8/3) hexagonal close-packed (hep) lattices 
p3| . The parameter a is the ratio of b to the dislocation core radius, r ; a ~ 2.9 for both 



bcc and fee crystals [|33[| . This melting relation plus experimental data on over half the 



elements in the periodic table give b 2 d(T m ) = 0.61 ±0.20 (throughout this paper the error 
in such expressions is the corresponding standard deviation) with G(300 °K), t;^s(300 °K) 
used instead of G(T m ), vws(T m ), respectively 



element 


Ba 


Cs 


Cr 


5-Fe 


K 


Li 


Na 


Nb 


Rb 


/3-Ti 


V 


T m , °K 


1000 


301.6 


2130 


1811 


336.5 


453.7 


370.9 


2750 


312.5 


1941 


2183 


v(T m ), A 3 


66.68 


116.8 


13.10 


12.76 


76.38 


22.14 


40.17 


19.33 


93.37 


18.61 


14.89 


G(T m ), GPa 


2.96 


0.39 


35.7 


30.8 


0.80 


3.60 


1.93 


32.6 


0.60 


21.9 


32.3 


Gv/(k B T m ) 


14.3 


10.9 


15.9 


15.7 


13.2 


12.7 


15.1 


16.6 


13.0 


15.2 


15.9 



Table 1. Numerical values of the ratio G(T m )v (T m )/(/c^T m ) for 11 elemental solids that 
melt from bcc crystalline structure at normal pressure. 



element 


Ag 


Al 


Ar 


An 


(3-Co 


Cu 


Ni 


Pb 


Pd 


Pt 


Rh 


T m , °K 


1235 


933.5 


83.8 


1338 


1768 


1358 


1728 


600.6 


1828 


2041 


2237 


v(T m ), A 3 


18.19 


17.55 


40.90 


17.88 


11.96 


12.61 


11.85 


31.14 


15.65 


16.04 


14.87 


G(T m ), GPa 


17.2 


15.6 


0.60 


15.2 


34.7 


27.1 


38.6 


5.60 


35.0 


32.0 


55.0 


Gv/(k B T m ) 


18.4 


21.2 


21.2 


14.7 


17.0 


18.2 


19.2 


21.0 


21.7 


18.2 


26.5 



Table 2. Numerical values of the ratio G(T m )v(T m )/(kBT m ) for 11 elemental solids that 
melt from fee crystalline structure at normal pressure. 



From the compilation of data in Tables 1 and 2 we find that the product of A and the 
logarithm in Eq. (2) (with v(T m ) = 0.42 ±0.02 p6| ) is a constant to 15% at zero pressure: 



A / o? 
8vr n \4b 2 d(T„ 



0.100 ±0.015, bcc, 
0.091 ±0.014, fee. 



(3) 



We make the reasonable assumption that the mean interdislocation distance at the 
melting point, 2R wl/ Jd(T m ), scales with b, which implies that b 2 d(T m ) is a compression- 
independent constant. It is also assumed that a' 1 = r /b is unchanged under compression. 
Hence A ln(a 2 /4& 2 <f) is expected to be pressure-independent with approximately the same 
value for both bcc and fee elements. It then follows from (2) that for a given element 



l-u(P,T m (P))/2 G(P,T m (P))v(P,T m (P)) 
1 - u(P, T m {P)) k B T m (P) hi(z - 1) 



(4) 



where the constant c has nearly the same value for both bcc and fee elements. Experi- 
mental validation of this relation is not posible because of a lack of data from moderate to 
high compressions. However, the P — > and P —>■ oo limits are consistent with Eq. (4), 
which we now demonstrate. 

At very high compressions a solid becomes a crystallized one-component plasma (OCP), 
i.e., a lattice of ions in a uniform neutralizing background of electrons [pSfl . The melting 
curve of a solid at ultrahigh pressures is described by the equation 



Z 2 e 2 



a(T m )k B T n 



(5) 



where Z is the atomic number, a = (3t>/47r) 1//3 is the Wigner-Seitz radius, and T m , a di 



mensionless constant, is the OCP coupling parameter at melt |3~5|. Numerous calculations 



of T m for a bcc OCP crystal (see ref. |37] for a review) converge on the value 175 |39f . 
The value of T m for a fee OCP crystal has been calculated to be 196 ± 1 |4TJ and 208.3 



41] ; hence we take T m = 200 for a fee OCP crystal in the following analysis. The bcc 
OCP single-crystal elastic constants (en — ci2)/2 and C44 have been calculated by means 
of Monte-Carlo simulations 42]| . A linear fit to the values of G given by the formula of 
Sisodia et al. |^3| (when c\\ and cu are not known separately, the value of G given by 
this formula approximates Kroner's shear modulus with high accuracy and, in fact, tends 
to the precise Kroner value in the limit P — > 00) results in |37] 



Cbcc P (T) - 9bcc ^— J 



MtA 1 / 3 Z 2 e 2 



A/3 



1 Mbcc rp 

J- m. 



(6) 



where sf bcc = 
coefficient g{ c 



0.09301 and /3^ p = 0.21 ± 0.18. We have calculated (unpublished) the 
to be 0.09011. The coefficient /?fc<F p has not been calculated, so we assume 



AccP P = Pbcc P ■ We have also calculated the Voigt (V) and Reuss (R) bounds on the shear 
modulus of an ideal hep OCP crystal: g^cp = 0.1194, g^ cp = 0.1045, hence ghc P = 0.1120 
for the Voigt-Reuss-Hill average. 

From Eqs. (4) and (6), and the ultrahigh pressure limit v{T) = 1/2 {44], 4~5| , we obtain 
focp = 9 9 ± 2.3 an d ^ocp = 8 9 ± 2 .o. Comparison of the OCP values of f to their 



zero-pressure counterparts (which follow from Eqs. (2) and (3)), £bcc(0) = 10.0 ± 1.5 and 
£fcc(0) = 11.0 ± 1.7, shows that the P = and OCP values agree to within uncertainties, 
compelling evidence, though not a proof, that Eq. (4) is in fact valid, at least for bcc and 
fee lattices. The uncertainty-weighted average of the four values is 10.0 ± 1.8. 

Formula (1) now follows from Eq. (4) provided that the ratio (l—u(T m )/2)/(l—v(T m )) 
is (approximately) a constant; in fact this ratio varies between w 4/3 at P = and 3/2 
as P -> oo, i.e., it is (17 ± 1)/12 « 17/12 to 94% accuracy. 

Formula (1) can also be derived from the theory of a Debye solid. Ledbetter |46[ 
derived the Debye-solid relation 

Q D = ^J~, (7) 
where Qd is the Debye temperature and A is a constant. (Since G ~ p 4//3 as p — > oo, 



counterpart 



which is consistent with 7 (Griineisen 



1/2 



01 



46 



Its widely used 

where B is the bulk modulus, has the wrong 



asymptotic behavior, 0£> ~ P ■) Siethoff and Ahlborn demonstrated the validity 
of the Ledbetter formula at P = for Debye-like cubic solids f|7|, |48|, |4"9|| , non-Debye 
hexagonal and tetragonal solids |50| , and intermetallic compounds fl51~| . Eq. (7), f ~ 1/p, 
and the Lindemann melting criterion |35| 



T m {p)p 



2/3 



constant, 



(8) 



again yield relation (1) 



T m {P), °K 


v{P,T m (P)), A 3 


u t , m/s 


G{P,T m {P)), GPa 


Gv/(k B T m ) 


205.59 


35.698 


952.6 


1.686 


21.21 


190.90 


36.216 


909.7 


1.516 


20.84 


175.91 


36.785 


879.5 


1.395 


21.14 


162.80 


37.319 


843.0 


1.263 


20.98 


162.07 


37.350 


847.0 


1.274 


21.28 


148.19 


37.959 


800.0 


1.118 


20.75 


134.47 


38.601 


768.6 


1.015 


21.11 


123.16 


39.155 


736.0 


0.918 


21.15 


83.80 


40.900 




0.600 


21.22 



Table 3. Numerical values of the ratio G(P,T m (P))v(P,T m (P))/(kBT m (P)) for Ar along 
its solidus. The last row of the table contains P = values. 



2.2 Experimental verification 



Direct experimental validation of relation (1) over a restricted range of densities is possible 
for a single element, viz. argon. Ishizaki et al. |52| measured the transverse ultrasonic 
wave velocity, u t , in compressed argon along its solidus as a function of temperature. 



We calculate the shear modulus from the formula ut = JG/p, and v = V/Na from the 



measured argon melting curve |53|, V = V(T m ) } V being the molar volume. Our results 
for the values of Gf/(fcsT m ) are shown in Table 3. 

For the P > data we find G v / '(fcsT m ) = 21.06 ± 0.17, in agreement with its P = 
value (we get Gv/(k B T m ) = 21.08 ± 0.17 for all of the data). Thus, Gv/(k B T m ) for Ar 
deviates from a constant by less than 1%. 



Model of the shear modulus at all temperatures 
and densities 



Preston and Wallace [gOJ constructed a model of the temperature dependence of the 
shear modulus (0 < T < T m ) for arbitrary pressures. The T-dependence of G involves two 
characteristic temperatures, namely the Debye temperature and the melting temperature. 
The shear modulus is always monotonically decreasing with T, and is nonlinear for T ~ <->£> 
and linear from 0£> to T m for most elements. An accurate representation of G(T) at fixed 
density is achieved by ignoring the low-temperature non-linearity and approximating G(T) 
as a linear function of the reduced temperature T/T m with the correct value G(p, 0) at 



T = p| 



G(p,T) = G(p, 0)(l-/?^-^, (9) 

In general, the parameter (3 may be density dependent. A fit to shear- modulus data 
spanning temperatures from T = to T/T m ~ 0.4 at zero pressure gave (3q = 0.23 ± 0.08 
pOf . (For the 11 fee elements in Table 4 below /3q = 0.27 ± 0.10.) On the other hand, 



(3 OCP = 0.21 ± 0.18 (as discussed above), which equals /3 to within uncertainties, so we 
assume that j3 is independent of density. At p = p ref and T = T m (p re f), Eq. (9) reduces to 

a -, _ G{p rc t,T m {p vc t)) 

P GW,0) ' 1 ' 

The linear temperature dependence is suggested by available P = experimental data 
on G over the temperature range < T < T m [1-19]. This straight-line representation 
turns out to be quite accurate: the maximum deviation of the data from the corresponding 



fitted lines is ~ 5 % for 21 of the 22 metals analyzed in ||30|| . The exception is uranium, 
for which G(T) is nonlinear throughout the a phase at P = 0. As mentioned above, G(T) 
is is nonlinear below 6^, thus G(T) is nonlinear for low-melting-point solids from T = 
to T m . Despite the nonlinearity of G(T) in these cases, the model uncertainty is only of 
order 10%. 

At any given pressure, the introduction of the temperature dependence of the density, 
p = p(T), into Eq. (9) gives the temperature dependence of G at that pressure. In Fig. 
1 we compare G(p(T), T) for < T < T m at P = for Au and Cu to experimental data 
0, |i~0"| . The temperature dependence of the density was taken from ref. fl54| , and G(p, 0) 



and T m (p) are described by Eqs. (13) and (14) below with parameter values from Tables 
2 and 4. 



200 400 600 800 1000 1200 

T, K 



Fig. 1. The P = shear moduli of Cu and Au: Eq. (9) with p = p{T) from ref. |54 
and G(p(T),0) and T m (p(T)) described by Eqs. (13) and (14) with the parameters from 
Tables 2 and 4 vs. the experimental data on Cu |H| (smaller points) and Au @ (larger 
points). 



The Griineisen parameter was recently modeled as J34 



1 ry-y 72 

lip) = o + "TTs + T> 7i> 72, q = const, q > 1, (11) 

through consideration of its low- and ultrahigh-pressure limits. This analytic form for 7 
was obtained under the assumptions that (i) 7 — > 1/2 as p — >• 00, (ii) 7 is an analytic 
function of x = 1/p 1 ^ 3 , essentially the interatomic distance, and (hi) the coefficient of x 
in the Taylor-Maclaurin series expansion for 7 is non-zero. The third term on the right- 
hand- side of Eq. (11) represents the contribution of the quadratic and higher-order terms 
in the power series. The procedure for calculating the values of 71, 72, and q is discussed 
below. 

Eq. (11) and the Lindemann criterion 

<"»T.(p) _ 2 / _ IN 



d In p V 3 

provide a model for the density dependence of the melting temperature 

T m (p)=T m (p Te{ )(-^J exp 1 671 ( (^1/3 ~ ~TH ) + — 1 77^7 ( 13 ) 




The natural choice for the reference density is p m , the zero-pressure density at melt, which 



is known experimentally in most cases (see, e.g., [54]). 
Finally, Eqs. (1),(9),(10), and (13) result in 

G(p , o) = G(p „„ 0) (^-) 4,S exp [ 6J1 - -k) + ^ - i) } , (14) 

where p re f is most conveniently chosen to be either p m or p , the density at zero pressure 
and temperature. 

Eqs. (9), (13), and (14) constitute our analytic model for the shear modulus. It requires 
the determination of 7 parameters, namely p re f, G(p re f,0), T m (p ref ), 72, q, and ft. 
The values of 71, 7 2 and q are obtained by simultaneous solution of Eq. (11) with p = 



p (T = 300 °K) and p = p m , and Eq. (5) with Y m = 180 |34j] and T m (p) given by the 
high-density limit of Eq. (13). The value of 7(p m ) is obtained from the Kraut-Kennedy 
relation and low-pressure melting data. The remaining parameters are either zero- 
pressure experimental data themselves or can be determined from such data (for example, 
ft). In Table 4 we present the values of p re f (both p and p m ), G(po,0), 71, 72, q, and 
ft for all of the fee elements of Table 2. The values of G(p m , 0) can be calculated from 
the relation G(p m , 0) = G(p m ,T m )/(l — ft) with G(p m ,T m ) from Table 2, which also 
contains the values of T m (p m ). Since /3-Co exists only above T ps 700 °K at P = 0, 
its values of G(p ,0) and ft were obtained from the conditions G(p m ,T m ) = 34.7 and 
G{p{T = 710 °K) = 8.62, T = 710 °K) = 57.1 pi. 



element 


Pa, g/cc 


Pm, g/cC 


G(p ,0), GPa 


7x, (g/cc) 1/3 


72, 


5/cc) 9 


Q 


ft 


Ag 


10.63 


9.850 


33.5 


2.23 


9.63 


10 4 


4.8 


0.18 


Al 


2.730 


2.550 


29.3 


0.84 


45.4 




3.5 


0.22 


Ar 


1.771 


1.622 


1.46 


1.06 


6.42 




2.2 


0.23 


Au 


19.49 


18.29 


30.5 


3.21 


1.97 


10 12 


9.4 


0.18 


/3-Co 


8.910 


8.180 


73.2 


1.81 


6.28 


10 4 


5.5 


0.33 


Cu 


9.020 


8.370 


52.4 


1.87 


2.31 


10 4 


4.7 


0.25 


Ni 


8.970 


8.220 


93.6 


1.85 


5.60 


10 5 


6.5 


0.41 


Pb 


11.60 


11.05 


11.7 


3.09 


8.21 


10 8 


8.5 


0.36 


Pd 


12.13 


11.29 


50.3 


2.40 


3.34 


10 6 


6.6 


0.07 


Pt 


21.58 


20.19 


66.3 


3.21 


1.13 


10 11 


8.3 


0.27 


Rh 


12.49 


11.49 


158. 


2.16 


1.46 


10 7 


6.5 


0.42 



Table 4. Numerical values of the model parameters for 11 fee elements. The correspond- 
ing values of T m (p m ) and G(p m , T m (p m )) are provided in Table 2. 



In Figs. 2, 3, and 4 we compare the melting curves of Ar, Al and Cu as given by 
Eq. (13) with the corresponding parameters from Table 4 to experimental data. 
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Fig. 2. Melting curve of Ar: Eq. (13) with the Ar parameters from Table 4 vs. data. The 
smaller points are the experimental data of ref . [[53| , and the larger points are the results 



of calculations 57 . 
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Fig. 3. Melting curve of Al: Eq. (13) with the Al parameters from Table 4 vs. data. The 
smaller points are the experimental data of ref. |5H] , and the larger points are the results 
of calculations 159 . 
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Fig. 4. Melting curve of Cu: Eq. (13) with the Cu parameters from Table 4 vs. data. 



The smaller points are from a new SESAME melting curve table for Cu ||60fl . The larger 
points are the P = reference point at (in g/cc) p = 8.4, and the shock-melting points 
of ref. H at p = 10.2, 12.2 and 14.3, and of refs. M, at p = 14.0. 



Only five of the seven shear modulus parameters are independent because four appear 
in the model as two ratios, namely (3/T m (p) in Eq. (9) and G(p re f, 0)/(p re f) 4 / 3 in Eq. (14); 
hence the shear modulus is of the form 

G(p, T) = di p 4//3 exp < - 3 — > — a 4 pT, a\, a,2, 03, a 4 , q = const > 0. (15) 

I P q P 1/3 J 

As specific examples we provide the following formulas for the shear moduli of Ar, Al, 
Cu, and Au: 

G Ar (p,T) = 687.4 p 4 / 3 exp j-^-^j - 1.32-KT 3 pT, (16) 

{ 25 9 5 04 1 
~-^J- ~ 1-85- 10- 3 pT, (17) 

G Cu (p,T) = 841.2 p 4 / 3 exp |-^^-^| -7.96- 10" 4 pT, (18) 

f 4 1 9 • 1 11 1 9 26 1 
G AU (P, T) = 1022.0 p 4 / 3 exp - ^ ^ - ^ - 1.37 • 10~* pT. (19) 



i20r 




2 . 5 



3 3.5 4 

rho, g/cc 
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Fig. 5. The T = 300 shear modulus of Ar: Eq. (16) vs. older |M] (smaller points) and 
more recent ||| (larger points) experimental data. The experimental technique used to 
obtain the older data has been criticized |^5| . 
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Fig. 6. The T = shear modulus of Al: Eq. (17) vs. the electronic-structure calculations 
of ref. |66| . The small, medium, and large points represent the values of G in fee, hep, 
and bec phases of Al, respectively. 



In Fig. 5 we compare the T = 300 °K shear modulus of Ar as given by Eq. (16) to 
experimental data. The T = shear modulus of Al from Eq. (17) is compared to the 
results of electronic structure calculations in Fig. 6. The T = shear moduli of Cu 
and Au as given, respectively, by Eqs. (18) and (19) are compared to the results of the 
corresponding electronic-structure calculations in Figs. 7 and 8 in the next section. 



4 Comparison of model to electronic-structure re- 
sults for Cu and Au 

With the exception of Ar, experimental data are not available to test the model to megabar 
pressures. We can, however, test the T = version of the model by comparing it to the 
results of ab initio electronic-structure calculations of the shear modulus. 

Electronic structure calculations based on approximate density functional theories 
have proven to give good predictions for a variety of material properties. A study of the 
elastic constants of several elements and compounds |J7| covering a wide range of elastic 
properties, found errors with respect to experiment of generally less than 10% in the 
isotropic shear modulus. These results are obtained without empirical inputs. We expect 
such calculations to have similar accuracy under compression, thus providing a test of the 
new analytic model. 

For this reason, we have carried out electronic structure calculations to obtain the 
single crystal elastic constants C = (Cu — Cu) /2, C44, and B — | (Cn + 2C12) for the 
fee metals Cu and Au from normal to twice normal density. From these an average 
polycrystalline shear modulus is calculated and compared to the model. 

The method for the calculations was described by Soderlind et al. pBfl . To evaluate 
C, the lattice is deformed by the (volume conserving) transformation 



(1 + 6 
1 + 5 
V 1/(1 + 6) 2 



(20) 



The resulting energy change is 



AE/V = QC'6 2 + 0(6 S ). 



(21) 



Similarly, C44 is obtained by applying the (volume conserving) transformation 

(16 
5 1 
V 1/(1 - 6 2 ) 



(22) 



resulting in an energy change 



AE/V = 2C 44 5 2 + 0{5 4 ). 



(23) 



In our calculations we have evaluated the energy as a function of 6 at intervals of 0.01, 
up to 6 = 0.04. For the C case the energy is not an even function of 6, and so negative 



values of delta were used. The resulting E(S) were fit to 4 th degree polynomials and the 
quadratic coefficent was used to evaluate the elastic constant from Eq. (]2l|) or Eq. (p3|). 

The bulk modulus B is obtained from the volume-dependent energy of the undistorted 
crystal by 

d 2 E , s 

B = V— y (24) 

The energy was evaluated at volume intervals of 5% of the normal volume, from 20% 
expanded to 50% contracted. Derivatives were evaluated by fitting the equation of state 
of Rose et al. p9j to the energies and differentiating the function. It was found that a 



single curve of this type did not accuaretly fit both the high density points and the points 
near the minimum, so seperate overlapping fits were made for the 10 highest and lowest 
densities. 

The electronic structure calculations were based on the linearized augmented plane- 
wave (LAPW) code WIEN97 JT[J. The energy functional used was the generalized- 



gradient approximation as parameterized by Perdew, Burke, and Ernzerhof |7TJ. Some 
numerical parameters used in the calculations for Cu (Au) were, in atomic units: muffin 
tin radius r MT = 1.8 (2.0), plane wave cut-off r M T^max = 9.0, cut-off for expansion of 
density and potential g max — 16 (19); Brillouin zone integrals used special points corre- 
sponding to 16 3 (18 3 ) points in the full zone, with Gaussian smearing of the energies by 
20 mRy 

Our results on C, C44 and B for Cu and Au are shown in Tables 5 and 6, respectively. 
It is interesting to note the increasing anisotropy of Au under pressure. From Table 6 we 
see that, for Au, C does not increase nearly as rapidly as C44 with compression. This 
is connected with the fact that the energy difference between the fee and bec structures 



is small at all pressures W%. The distortion corresponding to C is along the Bain path 



connecting fee to bec, and it has been seen |68| that a small energy difference between 



these structures correlates with a small value of C . 

Let us now turn to the calculation of the shear moduli of Cu and Au. For a solid 



of cubic crystalline structure, as analysis by Kroner |28j shows, successively narrower 
bounds can be placed on the shear modulus as the degree of disorder in grain orientation 
increases. In the limit of perfect disorder, the value of the shear modulus is the single 
positive real root of a cubic equation with coefficients that depend on the single-crystal 
elastic constants C, C 44 , and B : 

3 9£ + 4C" 2 3(B + 4C')C U 3BC'C U 
x H x — x — = 0. (25) 



The values of the shear modulus calculated from Eq. (^5J) are shown in Tables 5 and 
6 along with C, C44 and B. 

As a by-product of our analysis, we obtain the interesting results that G(2p Q , 0) ~ 
10G(p ,0) for Cu, and G(2p ,0) ~ 20G(p ,0) for Au. 

In Figs. 7 and 8 we compare Eqs. (18) and (19) with T = 0, for Cu and Au, to the 
corresponding G entries in Tables 5 and 6. 



p, g/cc 


C, GPa 


C44, GPa 


B, GPa 


G, GPa 


8.850 


30.404 


77.639 


142.15 


53.912 


9.833 


41.863 


124.78 


235.73 


81.901 


11.06 


48.599 


167.42 


386.02 


104.66 


12.64 


83.020 


260.81 


652.48 


168.57 


14.75 


130.48 


445.26 


1151.8 


279.70 


17.70 


229.71 


800.23 


2118.4 


499.25 



Table 5. The single-crystal elastic constants and shear modulus of Cu as functions of 
density from the electronic-structure calculations described in the text. 



p, g/cc 


C, GPa 


C44 , GPci 


B, GPa 


G, GPa 


19.29 


16.445 


31.690 


201.20 


24.585 


21.43 


19.550 


77.940 


339.57 


46.764 


24.11 


33.890 


127.75 


568.39 


78.093 


27.56 


35.053 


255.22 


1029.0 


127.48 


32.15 


69.837 


479.27 


1918.2 


243.34 


38.58 


121.16 


912.15 


3753.2 


451.65 



Table 6. The single-crystal elastic constants and shear modulus of Au as functions of 
density from the electronic-structure calculations described in the text. 
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Fig. 7. The T = shear modulus of Cu: Eq. (18) vs. electronic-structure calculations 
(larger points, Table 5). The smaller points, obtained from first-principles calculations 
p6| . are shown for comparison. 
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Fig. 8. The T = shear modulus of Au: Eq. (19) vs. electronic-structure calculations 
(larger points, Table 6). The smaller points, obtained from first-principles calculations 
|[73|| , are shown for comparison. 
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Fig. 9. Comparison of the two models for the T = shear modulus: Eqs. (16) and (17) vs. 
the corresponding Guinan-Steinberg values for Ar (smaller points) and Al (larger points). 



Finally, it is interesting to compare our model at T = to the Guinan-Steinberg 
model mentioned in the Introduction. The equation of state, P = P(p), is needed to 
make this comparison. In Fig. 9 the models are compared to each other for Ar and Al. 
The corresponding equations of state are taken from ref. [ff4[ . For Al, G' Q = 1.8 comes from 
ref. J75J. For Ar, G' = 1.6 is obtained from the relation 7 = B /2 G' /G - 
7o = 7(Po) from Eq. (11) and B taken from ref. [[F4[]. The values of p an d G 
can be found in Table 4. 

It is seen that agreement between the two models is good at low densities, but it 
gradually deteriorates with increasing compression. The reason for this must be that the 
Guinan-Steinberg model generally provides only the correct functional form G ~ p 4//3 in 
the limit of infinite compression, not the precise numerical value of G in that limit, in 
contrast to our model which provides both. 



1/6 with 

-G(p o ,0) 



5 Concluding remarks 

We have constructed an analytic model of the shear modulus applicable at all densities 
greater than or of order ambient (G(p, 0) — ► as p — ► 0, as required, but the model 
may not be quantitatively correct for expanded states), and temperatures from to T m . 
All of the model parameters can be obtained from low-pressure experimental data. The 
model has the proper low-pressure and high-pressure limits, by construction, and to within 
uncertainties it agrees with electronic-structure values of G for Cu and Au to compressions 
of 2, which roughly corresponds to pressures of 5 Mbar for Cu and 9 Mbar for Au. 

The above comparisons of our shear modulus model, which includes a model for T m (p), 
to electronic structure calculations and experimental data on Ar, Al, Cu, and Au show 
very good agreement. This suggests that our model accurately represents the density 
and temperature dependence of the shear moduli of monatomic solids in general. There 
is, however, no theoretical justification for applying our model to alloys or compounds, 
although in practice it may work reasonably well in these cases. Its generalization to more 
complex materials would involve generalizing our model for the Gruneisen parameter. 
A functional form for ~f(p) depends critically on the asymptotic (p — > oo) form of the 
equation of state |33| , and it has been suggested that the asymptotic forms of the equations 
of state of more complex materials, e.g., ionic, covalent, or molecular crystals, are different 



from that of a metal|76|. If so, the limiting value of 7 is unknown (not necessarily 1/2) 
for such materials. In that case, an analytic model for the Gruneisen parameter cannot 
be constructed, hence analytic forms for the melting curve and shear modulus cannot be 
obtained. 

We now briefly discuss three potential applications of our model. 

(1) Plastic deformation of metals at high pressure. It is generally assumed that the 
ratio of the plastic flow stress (shear stress necessary to induce plastic deformation at 
a given strain rate) to the shear modulus is approximately independent of pressure. In 
other words, the predominant pressure dependence of the plastic flow stress is contained 
in the shear modulus. An accurate, simple analytic (for fast evaluation) model of the 
shear modulus is therefore essential for numerical simulations of material deformation 
over extremes in pressure. 



(2) Numerical simulations of elastic wave propagation, including pressure release waves 
in shocked solids. The differential stress deviator, ds^j, is equal to 2G(p,T) (de^ — 
$ij de^k/Si) plus material rotation terms (de^ is the differential elastic strain), thus a model 
of the shear modulus is required for calculations of elastic wave propagation in materials 
with sufficiently high yield stresses that the stress deviators are not negligible. The speed 
of a release wave in a shocked solid depends on G(ph,T h ), where pu and T H are the 
density and temperature of the shocked state. 



(3) Pulsations and quakes of dense stars. Hansen and Van Horn [77] have done a 
preliminary analysis of the effects of crystalline cores on the oscillations of white dwarfs 
and found that the g-like spheroidal mode frequencies are increased by approximately a 
factor of two, concluding that the elastic shear strength of the core must be taken into 
account in the computation of cool white dwarf oscillations. The inclusion of elastic shear 



strength in the neutron star pulsation equations of McDermott et al. |78| resulted in the 
appearance of two classes of oscillation modes not present in a fluid neutron star. The 
change in the shape of the surface following a neutron star quake is proportional to the 



shear modulus of the crust [79 



Further tests of our model for the shear modulus should be made as high-pressure 
experimental data and electronic structure results become available for elements other 
than argon, aluminum, copper and gold. New zero-pressure data are also needed to 
generate additional sets of model parameters. 
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